require(mvtnorm)
require(OpenMx)
require(psych)
require(MASS)

rAC  <- 0                                                # Set the level of the correlation betwween the A and C Varaince Components
nsim <- 200000                                           # Set the number of simulations (We used 200000)

MZdatasets <- list(NULL)
DZdatasets <- list(NULL)

ests <- matrix(NA, nsim, 6)


nv <- 1                                                     # Set the number of variables that you will be using (per person)
Vars <- c('t')                                              # Name the variable
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")        # Create an object that will select the variables
ntv <- nv*2                                                 # Set the total number of variable that you will use (1 var * 2 twins)




mzmat <- matrix(c(  1,  1,rAC,rAC,  0,  0,
                    1,  1,rAC,rAC,  0,  0,
                  rAC,rAC,  1,  1,  0,  0,
                  rAC,rAC,  1,  1,  0,  0,
                    0,  0,  0,  0,  1,  0,
                    0,  0,  0,  0,  0,  1),6)                       # Create a matrix of correlations for the genetic and environmental variance components for each twin (MZ)


dzmat <- matrix(c(  1,1/2,rAC,rAC,  0,  0,
                  1/2,  1,rAC,rAC,  0,  0,
                  rAC,rAC,  1,  1,  0,  0,
                  rAC,rAC,  1,  1,  0,  0,
                    0,  0,  0,  0,  1,  0,
                    0,  0,  0,  0,  0,  1),6)                       # Create a matrix of correlations for the genetic and environmental variance components for each twin (DZ)


 for(sim in 1:nsim) {

	 mzSim  <- mvrnorm(n=1000, c(0,0,0,0,0,0), mzmat , empirical=T)   # Simulate genetic and environmental factor scores for each MZ twin
	 dzSim  <- mvrnorm(n=1000, c(0,0,0,0,0,0), dzmat , empirical=T)   # Simulate genetic and environmental factor scores for each DZ twin

VarCompE <- runif(1,.2,.5)                                            # Set the Genetic Variance Component
VarCompA <- runif(1,0,(1-VarCompE))                                   # Set the Common Environmental Variance Component
VarCompC <- 1-VarCompA -VarCompE                                      # Set the Unique Environmental Variance Component

	 aPath <- sqrt(VarCompA)                                          # Set the Genetic Variance Path Coefficeint
	 cPath <- sqrt(VarCompC)                                          # Set the Common Environmental Path Coefficeint
	 ePath <- sqrt(VarCompE)                                          # Set the Unique Environmental Path Coefficeint

	 MZdata <- cbind((mzSim[,1]*aPath+
	                  mzSim[,3]*cPath+
	                  mzSim[,5]*ePath),
	                 (mzSim[,2]*aPath+
	                  mzSim[,4]*cPath+
	                  mzSim[,6]*ePath))                                # Construct the Phenotypes based on genetic and environmental factor scores MZ twins

	 DZdata <- cbind((dzSim[,1]*aPath+
	                  dzSim[,3]*cPath+
	                  dzSim[,5]*ePath),
	                 (dzSim[,2]*aPath+
	                  dzSim[,4]*cPath+
	                  dzSim[,6]*ePath))                                # Construct the Phenotypes based on genetic and environmental factor scores for DZ twins

	 colnames (MZdata) <- selVars
	 colnames (DZdata) <- selVars
		
		MZdatasets[[sim]] <- MZdata		                                           # Save the MZ data
		DZdatasets[[sim]] <- DZdata		                                           # Save the DZ data
		ests[sim,1:6] <- c(aPath,cPath,ePath,VarCompA,VarCompC,VarCompE)       # Save the simulated parameter values for each rep
		
	}
	colnames(ests) <- c("aPath","cPath","ePath","VarCompA","VarCompC","VarCompE")   

#MZdatasets
#DZdatasets